Stochastic modeling
of reaction networks is a framework used to describe the time evolution
of many natural and artificial systems, including, biochemical reactive
systems at the molecular level, viral kinetics, the spread of epidemic
diseases, and wireless communication networks, among many other
examples. In this work, we present a novel multilevel Monte Carlo method
for kinetic simulation of stochastic reaction networks that is speci
cally designed for systems in which the set of reaction channels can be
adaptively partitioned into two subsets characterized by either "high"
or "low" activity. Adaptive in this context means that the partition
evolves in time according to the states visited by the stochastic paths
of the system. To estimate expected values of observables of the system
at a prescribed final time, our method bounds the global computational
error to be below a prescribed tolerance, TOL, within a given confidence
level. This is achieved with a computational complexity of order O (TOL^{-2})
, the same as with an exact method, but with a smaller constant. We
also present a novel control variate technique based on the stochastic
time change representation by Kurtz, which may dramatically reduce the
variance of the coarsest level at a negligible computational cost. Our
numerical examples show substantial gains with respect to the standard
Stochastic Simulation Algorithm (SSA) by Gillespie and also our previous
hybrid Chernoff tau-leap method.